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Abstract: We investigate the radiation pressure feedback in the formation of massive stars in 1, 2, and 3D radi- 
ation hydrodynamics simulations of the collapse of massive pre-stellar cores. In contrast to previous research, 
we consider frequency dependent stellar radiation feedback, resolve the dust sublimation front in the vicinity 
of the forming star down to 1.27 AU, compute the evolution for several 10 5 yrs covering the whole accretion 
phase of the forming star, and perform a comprehensive survey of the parameter space. The most fundamental 
result is that the formation of a massive accretion disk in slowly rotating cores preserves a high anisotropy in 
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the radiation field. The thermal radiation escapes through the optically thin atmosphere, effectively diminishing 
the radiation pressure feedback onto the accretion flow. Gravitational torques in the self-gravitating disk drive 
i— —i a sufficiently high accretion rate to overcome the residual radiation pressure. Simultaneously, the radiation 
pressure launches an outflow in the bipolar direction, which grows in angle with time and releases a substantial 
J> fraction of the initial core mass from the star-disk system. Summarized, for an initial core mass of 60, 120, 240, 
and 480 M these mechanisms allow the star to grow up to 28.2, 56.5, 92.6, and at least 137.2 M respectively. 
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The understanding of massive stars suffers from the lack of a generally accepted formation scenario. 
If the formation of high-mass stars is treated as a scaled-up version of low-mass star formation, a 
special feature of these high-mass proto-stars is the interaction of the accretion flow with the strong 
irradiation emitted by the newborn stars due to their short Kelvin-Helmholtz contraction timescale 
(Shu et al. 1987). Early one-dimensional studies (e.g. Larson & Starrfield 1971; Kahn 1974; Yorke 
) agree on the fact that the growing radiation pressure potentially stops and reverts the 
accretion flow onto a massive star. 

But this radiative impact strongly depends on the geometry of the stellar environment (Nakano 
59). The possibility was suggested to overcome this radiation pressure barrier via the formation of 
a long-living massive circumstellar disk, which forces the generation of a strong anisotropic feature 
of the thermal radiation field. Earlier investigations by Yorke & Sonnhalter (2002) tried to identify 
such an anisotropy, but their simulations show an early end of the disk accretion phase shortly after 
its formation due to strong radiation pressure feedback. Krumholz et al. (2( )9) demonstrated the 
possibility to maintain the shielding property of the accretion disk by additional feeding due to the 
so-called '3D radiative Rayleigh-Taylor instability' in the outflow cavity. 
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2 Method 



2.1 The aim 

Aim of this study is to reveal the details of the radiation dust interaction in the formation of the most 
massive star during the collapse of a pre-stellar core (gravitationally unstable fragment of a molecular 
cloud). Due to the strong focus onto the core center, the physics in the outer core region have to be 
simplified, e.g. further fragmentation of the pre-stellar core is suppressed. 



2.2 The code 

For this purpose, we use our newly developed self-gravity radiation hydrodynamics code. The hydro- 
dynamical evolution is computed using the open source code Pluto3 (Mignone et al. 2007), including 
full tensor viscosity. The frequency dependent hybrid radiation transport method is summarized by 
Kuiper et al. (2010a). The implementation of Poisson's equation, the viscosity prescription as well as 
the tabulated dust and stellar evolution model are presented in Kuiper et al. (2010b). 

The simulations are performed on a time independent grid in spherical coordinates. The radially 
inner and outer boundary of the computational domain are semi-permeable walls, i.e. the gas can 
leave but not enter the computational domain. The resolution of the non-uniform grid is chosen to be 
(Ar x r A6 x r sin(0) A0) min = 1.27 AU x 1.04 AU x 1.04 AU in the midplane (6 = 90°) around 
the forming massive star and decreases logarithmically in the radial outward direction. The outer core 
radius is fixed to r max = 0.1 pc. The accurate size of the inner sink cell was determined in a parameter 
scan presented in Kuiper et al. (2010b), Sect. 5.1. 



2.3 The initial conditions 

The basic physical initial conditions are very similar to the one used by Yorke & Sonnhalter (2002). 
We start from a cold (T = 20 K) pre-stellar core of gas and dust. The initial dust to gas mass 
ratio is chosen to be M dust /M gas = 1%. The model describes a so-called quiescent collapse sce- 
nario without turbulent motion (u r — uq = 0). The core is initially in slow solid-body rotation 
(lu^l/R = Vl = 5 x 10~ 13 Hz). The initial density slope drops with r~ 2 and the total mass M core 
varies in the simulations from 60 up to 480 M© . 
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Figure 1 : Accretion rate M* as a function of the 
stellar mass M* for varying initial core masses. 
The periods of accretion are mentioned for each 
run. The vertical line marks the upper mass limit 
for spherically symmetric accretion flows. 



Figure 2: Accretion history of a three-dimensional 
run, driven by evolving torques in the self- 
gravitating massive accretion disk, compared with 
the axially symmetric viscous accretion models 
with and without a-viscosity (M core = 120 M ). 



3 Results 



The presentation and discussion of the selective simulation results are divided in two subsections, 
which distinguish between the radiation pressure feedback onto the disk region and its interaction 
with the accretion flow (Sect. 3.1), and the feedback onto the region of the bipolar cavities yielding 
the launching of stable radiation pressure driven outflows (Sect. 3.2). 

3.1 Radiation pressure in the circumstellar accretion disk 

First, we performed pre-stellar core collapse simulations in axial symmetry for four different initial 
core masses M core = 60 M , 120 M , 240 M , and 480 M . The resulting accretion histories as 
a function of the actual stellar mass are displayed in Fig. 1. Contrary to one-dimensional in-fall 
simulations (e.g. Larson & Starrfield 1971; Kahn 1974; Yorke & Kriigel 1977; Kuiper et al. 2010b), 
the disk accretion scenario allows the formation of massive stars far beyond the upper mass limit of 
spherically symmetric accretion flows (the so-called 'radiation pressure barrier'). The reason for this 
strong dependence of the radiation pressure feedback on the geometry of the stellar environment is 
that the massive accretion disk results in a high anisotropy of the radiation field. Most of the radiation 
escapes perpendicular to the accretion flow. Consistently, all simulations performed so far (including 
resolution studies and different sizes of the inner sink cell) show that the accretion rate increases with 
a stronger suppression of the isotropic part of the thermal radiation field. 

The difference between our simulations and the ones performed by Yorke & Sonnhalter (2002), 
which yield an early and abrupt end of the disk accretion phase directly after the formation of the 
circumstellar disk, is the fact that in our simulations the radiation dust interaction is resolved down to 
the dust sublimation front of the forming massive star. E.g. in contrast to the size of the inner sink cell 
of ?"mm = 80 AU or 160 AU as well as a resolution of the inner core region of (Ar) min = 80 AU or 
160 AU in the 60 M or 120 M collapse by Yorke & Sonnhalter (2002) respectively, we use an inner 
rim of the computational domain of r min = 10 AU with a maximum resolution of (Ar) min = 1.27 AU. 
For comparison: Radiation pressure becomes as strong as gravity for roughly a 20 M star (the so- 
called 'Generalized Eddington barrier') and its dust sublimation radius in the disk's midplane is about 
20 to 30 AU. The short disk accretion phases by Yorke & Sonnhalter (2002) are reproduced, if we use 
the same huge inner sink cell radii they used (see Kuiper et al. 2010b). 

Secondly, the validity of the standard a-shear-viscosity prescription based on Shakura & Sunyaev 
(1973), which we use in the axially symmetric disk simulations, is cross-checked in three-dimensional 
simulations without shear-viscosity, in which the angular momentum transport is self-consistently 
driven by evolving gravitational torques in the self-gravitating massive accretion disks. A comparison 
of the stellar mass growth after the onset of disk formation is shown in Fig. 2 for the collapse of a 
120 M pre-stellar core. Although the three-dimensional run yields a more time dependent episodic 
accretion history than the viscous disk simulations, the mean accretion rates equal each other. We 
conclude that the self-gravity of massive accretion disks is a sufficient driver of angular momentum 
transport. The resulting accretion flow easily overcomes the remnant thermal radiation pressure, 
which is - deeply inside the disk - diminished by the shielding effect at the inner rim of the dust disk. 

The difference between our simulations and Krumholz et al. (2007, 2009) is the fact that in the disk 
accretion scenario, which we propose here, the most massive stars can form consistently with standard 
accretion disk theories without the requirement of '3D radiative Rayleigh-Taylor instabilities' in the 
outflow region. We discuss the issue of the outflow stability in more detail in the following Sect. 3.2. 



3.2 Radiation pressure in the bipolar outflow 

The feedback of the radiation pressure onto the bipolar region differs strongly from the feedback onto 
the circumstellar disk region. The strong radiation pressure in the bipolar direction stops the in-fall 
and reverts the mass flow into a radiation pressure driven outflow with velocities of a few 100 km s _1 . 
Fig. 3 visualizes a slice through the three-dimensional density structure of such an outflow. The 




Figure 3: A slice through the three-dimensional density structure of a radiation pressure driven out- 
flow during a 120 M pre-stellar core collapse. 

radiation pressure driven outflows in our axially symmetric simulations remain stable throughout 
the whole stellar accretion phase. As a result, the radiation pressure in the bipolar direction thrusts 
a substantial fraction (more than 40% during the stellar disk accretion phase) of the initial mass 
of the pre-stellar core into outer space, reducing the star formation efficiency of the massive star 
forming core. In contrast to Krumholz et al. (2009), the outflows in our three-dimensional simulations 
remain stable so far as well in spite of the strong non-axially symmetric features developing in the 
outflow region. The reason(s) for that disagreement can be manifold, and we briefly overview the 
most relevant explanations here: 

1. It was supposed by Krumholz et al. (2009) that the instability requires non-axially symmetric 
modes, and this could explain why it is not seen in our two-dimensional simulations. In our 
three-dimensional runs, the outflow is stable during the launching and the first growth phase, but 
instability might appear at a later epoch. The instability in Krumholz et al. (2009) occurs when the 
outflow extends up to r ~ 1000 AU, and the outflow in our runs have propagated up to r < 500 AU 
so far. Our three-dimensional simulations are continuously running, and we will further report on 
their status in the near future. 

2. Our simulations account for frequency dependent stellar irradiation instead of the frequency aver- 
aged (gray) approximation used tl. (2009). 

3. We use a more accurate treatment (ray-tracing) of the first absorption of the stellar irradiation and 
the instability could be an artifact of the flux-limited diffusion approximation, which results in a 
wrong radiative flux at t ~ 1 directly in the cavity wall. 

4. The disagreement could be due to the different numerical resolution of the diverse grid structures. 

Conclusion: The details of the underlying physics of the '3D radiative Rayleigh-Taylor instability' 
are quite unknown so far and deserve further investigation. 



4 Summary 



We performed high-resolution radiation hydrodynamics simulations of monolithic pre-stellar core 
collapses including frequency dependent radiative feedback. The dust sublimation front of the form- 
ing star could be resolved down to 1.27 AU. The frequency dependent ray-tracing of our newly de- 
veloped hybrid radiation scheme denotes the most realistic radiation transport method used in multi- 
dimensional hydrodynamics simulations of massive star formation by now. In axial symmetry, the 
whole accretion phase of several 10 5 yrs was computed. The comprehensive parameter studies reveal 
new insights of the radiative feedback onto the accretion flow during the formation of a massive star: 

The massive accretion disk induces a strong anisotropy of the thermal radiation field. In contrast 
to the short disk accretion phases in Yorke & Sonnhalter (2002), the accretion disks reveal a persistent 
anisotropy also during their long-term evolution due to the fact that we resolve the dust disk down to 
the dust sublimation front of the forming star. Contrary to this explanation, Krumholz et al. (2009) 
stated that the much longer accretion phases in their own simulations compared to Yorke & Sonnhalter 
(2002) are a result of the so-called '3D radiative Rayleigh-Taylor instability' in the outflow cavity; 
in axially symmetric simulations sufficient feeding of the disk should therefore be impossible. But 
our axially symmetric collapse simulations show a stable radiation pressure driven outflow, while the 
circumstellar disk gains enough mass from the in-falling envelope to maintain its shielding property 
over several free fall times, in fact over a longer period ever simulated in previous research studies. 

Our three-dimensional core collapse simulations support this accretion disk scenario consistently. 
Gravitational torques in the self-gravitating disk lead to a sufficiently high angular momentum trans- 
port to maintain accretion. The accretion history is stronger episodically compared to the viscous disk 
evolution in the axially symmetric runs, but adds up to a very similar mean accretion rate. 

In the bipolar direction, the radiative feedback yields the launch of a stable radiation pressure 
driven outflow. No '3D radiative Rayleigh-Taylor instability' is identified yet. Although the three- 
dimensional simulations are still running, we can certainly conclude from the axially symmetric runs 
that additional feeding by an unstable outflow is not necessary to allow the formation of massive stars. 
On the other hand, a well-grounded conclusion on the stability of the outflow and its influence on the 
formation of massive stars clearly requires further investigation. 

Summing up, the central stars in our simulations of the disk accretion scenario grow far beyond 
the upper mass limit found in the case of spherically symmetric accretion flows. For an initial mass 
of the pre-stellar host core of 60, 120, 240, and 480 M the masses of the final stars formed add up 
to 28.2, 56.5, 92.6, and at least 137.2 M Q respectively. Indeed, the final massive stars are the most 
massive stars ever formed in a multi-dimensional radiation hydrodynamics simulation so far. 
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